% script to test Ute's data screening methods
% latest date 25.11.2004


function[] = ute(SPMatrix,ProcDef),
global Exx_p Ey_p Hx_p Hy_p Z ExPC EyPC Epol Hpol grey;
grey = [0.6 0.6 0.6];
[tf2,ProcDef,SPM2] = robspm(SPMatrix,ProcDef);
flag = 1;
i = 1;
nlen = size(SPMatrix,2);
LL = ones([size(SPMatrix,1),nlen]) == 1; % stupid way of getting boolean 1 ones
%L(300:400) = 0;
%global Ex_p,Ey_p,Hx_p,Hy_p;
xdata = 1:nlen;
while flag == 1,
    L = LL(i,:);
    [Exx_p,Ey_p,Hx_p,Hy_p]=getpower(SPMatrix,i);% PSD of fields
    Z = tf_dn_new(SPMatrix,ProcDef,i); % transfer functions & variance [Zxy,Zxy,Zxx,Zyy,ZxyVar....]
    [Epol,Hpol] = polangle(SPMatrix,i); % polarization angles
    Epol = Epol .* (180/pi); 
    Hpol = Hpol .* (180/pi);
    
    ExPC = pcoh(squeeze(SPMatrix(i,:,:,:)),'Ex');
    EyPC = pcoh(squeeze(SPMatrix(i,:,:,:)),'Ey');
    
    subplot(3,3,1);
    semilogy(Exx_p,'r.');
    axis([1 nlen -inf inf]);
    if sum(~L) > 0,
        hold on;
        semilogy(xdata(~L),Exx_p(~L),'m.');
        hold off;
    end;
    title('Power of Ex');
    subplot(3,3,2);
    semilogy(Hx_p,'r.');
     if sum(~L) > 0,
        hold on;
        semilogy(xdata(~L),Hx_p(~L),'m.');
        hold off;
    end;
    title('Power of Hx');
    subplot(3,3,3);
    semilogy(Hy_p,'r.');
    if sum(~L) > 0,
        hold on;
        semilogy(xdata(~L),Hy_p(~L),'m.');
        hold off;
    end;
    title('Power of Hy');
    subplot(3,3,4);
    plot(real(Z(:,1)),imag(Z(:,1)),'r.');
    hold on;
    plot(real(Z(:,3)),imag(Z(:,3)),'b.');
    plot(real(mean(Z(:,3))),imag(mean(Z(:,3))),'w*');
    plot(real(mean(Z(:,1))),imag(mean(Z(:,1))),'w*');
     if sum(~L) > 0,
        hold on;
        plot(real(Z(~L,1)),imag(Z(~L,1)),'m.');
    end;
    title('Impedance Re Vs Im');
     if sum(L) > 2,
        axis([min(real(Z(L,1))) max(real(Z(L,1))) min(imag(Z(L,1))) max(imag(Z(L,1)))]);
    else,
        axis equal;
    end;
    hold off;
    subplot(3,3,5);
    plot(Z(:,5),'r.');
    hold on;
    plot(Z(:,8),'b.');
    if sum(~L) > 0,
        hold on;
        plot(xdata(~L),Z(~L,5),'m.');
        plot(xdata(~L),Z(~L,8),'m.');
    end;
    legend('ZxyErr','ZxxErr');
    title('Bivariance  Error of Zxy & Zxx');
    hold off;
    subplot(3,3,6);
    plot(ExPC,'r.');
    if sum(~L) > 0,
        hold on;
        plot(xdata(~L),ExPC(~L),'m.');
    end;
    hold off;
    title('Predicted Coherence Ex to HxHy');
    subplot(3,3,8);
    plot(Epol,'r.');
    hold on;
    plot(Hpol,'b.');
    if sum(~L) > 0,
        plot(xdata(~L),Epol(~L),'m.');
        plot(xdata(~L),Hpol(~L),'m.');
    end;
    title('Polarization Angles');
    legend('E-pol','H-pol');
    hold off;
    disp1(tf2,ProcDef);
%   tf = tf_dn(squeeze(SPMatrix(i,:,:,:)));
    disp2(Z,ProcDef.TLFreq1(i),L);
%     h = axes('Position',[0 0 1 1],'Visible','off');
%     str(1) = {ProcDef.file1};
%     str(2) = {ProcDef.nband};
%     str(3) = {ProcDef.TLFreq1(i)};
%     text(0.1,0.25,str);
    fprintf('Current sum L = %10d\n', sum(LL(i,:)));

    men1 = menu('SEL','->','Exit','Sel','<-','DeSel','AllSel','AllDes','RunRob');
    if men1 == 1 & i < length(ProcDef.TLFreq1),
        i = i + 1;
    elseif men1 == 2,
        flag = 0;
    elseif men1 == 4 & i > 1,
        i = i - 1;
    elseif men1 == 3,
        LL(i,:) = updateL(L,1);
        fprintf('Updated sum L = %10d\n', sum(LL(i,:)));
    elseif men1 ==5,
        LL(i,:)  = updateL(L,0);
        fprintf('Updated sum L = %10d\n', sum(LL(i,:)));
    elseif men1 == 6,
        LL(i,:) = 1;
    elseif men1 == 7,
        LL(i,:) = 0;
    elseif men1 == 8,
        ProcDef.len = sum(LL(i,:));
        [tf2,ProcDef,SPM2] = robspm(SPMatrix(:,LL(i,:),:,:),ProcDef);
    end;
        
       
end;

function[Ex,Ey,Hx,Hy]=getpower(SPMatrix,nfrq),

for i = 1:length(SPMatrix),
    Hx(i) = SPMatrix(nfrq,i,1,1);
    Hy(i) = SPMatrix(nfrq,i,2,2);
    Ex(i) = SPMatrix(nfrq,i,4,4);
    Ey(i) = SPMatrix(nfrq,i,5,5);
end;



function[] = disp1(tf,ProcDef),

frq=ProcDef.TLFreq1;
file=ProcDef.file1;
nfrq = length(frq);

for i = 1:nfrq,   
   Rhxy(i) = 0.2 * (1/frq(i)) * (tf(i,1) * conj(tf(i,1)));
   Rhyx(i) = 0.2 * (1/frq(i)) * (tf(i,2) * conj(tf(i,2)));
end;

RhxyErr = (tf(:,5).*abs(tf(:,1))./frq') * 0.4; % bit confusion. may be ok -dt 2.8.2
RhyxErr = (tf(:,6).*abs(tf(:,2))./frq') * 0.4;

%--------------------------------------------------------------
T = 1./frq;

mnx=ceil(log10(max(T)));
mux=floor(log10(min(T)));
mny=ceil(log10(max([Rhxy Rhyx])));
muy=floor(log10(min([Rhyx Rhxy])));

c1 = 'r*';
c2 = 'b*';
c3 = 'red';
c4 = 'blue';

subplot(3,3,9);
%plot(log10(T'),log10(Rhxy),'r.',log10(T),log10(Rhyx),'b.',log10(T'),log10(Rhxy),'r-',log10(T),log10(Rhyx),'b-');

if ishold == 1,
   c1 = 'c*';
   c2 = 'g*';
   c3 = 'cyan';
   c4 = 'green';
end;

  % plot(log10(T'),log10(Rhxy),c1,log10(T),log10(Rhyx),c2);
   plot(log10(T'),log10(Rhxy),c1); %only plots the rhxy 
for i = 1:nfrq,
   tt = log10(Rhxy(i));
   kk = log10(Rhxy(i)+RhxyErr(i)/2);
   vv = 2*tt-kk;
   h=line(log10([T(i) T(i)]),[kk vv]);
   set(h,'color',c3);
%    tt = log10(Rhyx(i));
%    kk = log10(Rhyx(i)+RhyxErr(i)/2);
%    vv = 2*tt-kk;
% 
%   	h=line(log10([T(i) T(i)]),[kk vv]);
%    set(h,'color',c4);

end;
axis([mux mnx muy mny]);
grid;
ylabel('log ohm-m');
xlabel('log period (s)');
legend('XY','YX');
hold on;

function[] = disp2(tf,frq_c,L),
global Rhxy Rhyx grey;
Rhxy = [];
Rhyx = [];
nlen = length(tf);

for i = 1:nlen,   
   Phxy(i) = atan2(imag(tf(i,1)),real(tf(i,1))) * 180 / 3.1415926535897;
   Rhxy(i) = 0.2 * (1/frq_c) * (tf(i,1) * conj(tf(i,1)));
   Phyx(i) = (atan2(imag(tf(i,2)),real(tf(i,2))) * 180 / 3.1415926535897);
   Rhyx(i) = 0.2 * (1/frq_c) * (tf(i,2) * conj(tf(i,2)));
end;

T = 1./(ones([1,nlen])*frq_c);

c1 = 'r*';
c2 = 'b*';

subplot(3,3,9);
%plot(log10(T'),log10(Rhxy),'r.',log10(T),log10(Rhyx),'b.',log10(T'),log10(Rhxy),'r-',log10(T),log10(Rhyx),'b-');
%plot(log10(T'),log10(Rhxy),c1,log10(T),log10(Rhyx),c2); 
plot(log10(T'),log10(Rhxy),c1);
hold on;
plot(log10(T(~L)'),log10(Rhxy(~L)),c2);
hold off;
%-----------------------------------------------------------------------
function[L] = updateL(L,flag),
global Exx_p Ey_p Hx_p Hy_p Z ExPC EyPC Epol Hpol Rhxy;
men1 = menu('Sel Param','PowEx','PowHx','PowHy','Zxy','Biv Err','Pred Coh','PAngle','Rhxy');

if men1 == 1,
    subplot(3,3,1);
    [xmin,xmax,ymin,ymax] = getselection;
    L = TrimData(Exx_p,L,xmin,xmax,ymin,ymax,flag);
    fprintf('%10.0f %10.0f %10.3e %10.3e sum L = %10d\n', xmax,xmin,ymax,ymin,sum(L));
elseif men1 == 2,
    subplot(3,3,2);
     [xmin,xmax,ymin,ymax] = getselection;
      L = TrimData(Hx_p,L,xmin,xmax,ymin,ymax,flag);
    fprintf('%10.0f %10.0f %10.3e %10.3e sum L = %10d\n', xmax,xmin,ymax,ymin,sum(L));

  elseif men1 == 3,
    subplot(3,3,3);
    [xmin,xmax,ymin,ymax] = getselection;
    L = TrimData(Hy_p,L,xmin,xmax,ymin,ymax,flag);
    fprintf('%10.0f %10.0f %10.3e %10.3e sum L = %10d\n', xmax,xmin,ymax,ymin,sum(L));
  elseif men1 == 4,
    subplot(3,3,4);
    [xmin,xmax,ymin,ymax] = getselection;
    L = TrimData(Z(:,1)',L,xmin,xmax,ymin,ymax,flag);
    fprintf('%10.0f %10.0f %10.3e %10.3e sum L = %10d\n', xmax,xmin,ymax,ymin,sum(L));
  elseif men1 == 5,
    subplot(3,3,5);
    [xmin,xmax,ymin,ymax] = getselection;
    L = TrimData(Z(:,5)',L,xmin,xmax,ymin,ymax,flag);
    fprintf('%10.0f %10.0f %10.3e %10.3e sum L = %10d\n', xmax,xmin,ymax,ymin,sum(L));
  elseif men1 == 6,
    subplot(3,3,6);
    [xmin,xmax,ymin,ymax] = getselection;
    L = TrimData(ExPC,L,xmin,xmax,ymin,ymax,flag);
    fprintf('%10.0f %10.0f %10.3e %10.3e sum L = %10d\n', xmax,xmin,ymax,ymin,sum(L));
  elseif men1 ==7,
    subplot(3,3,8);
    [xmin,xmax,ymin,ymax] = getselection;
    L = TrimData(Hpol,L,xmin,xmax,ymin,ymax,flag);
    fprintf('%10.0f %10.0f %10.3e %10.3e sum L = %10d\n', xmax,xmin,ymax,ymin,sum(L));
elseif men1 == 8,
    subplot(3,3,9);
    [xmin,xmax,ymin,ymax] = getselection;
    L = TrimData(Rhxy,L,1,length(Rhxy),10^ymin,10^ymax,flag);
    
  end;
%-----------------------------------------------------------------------
function[xmin,xmax,ymin,ymax] = getselection,
    k = waitforbuttonpress;
    point1 = get(gca,'CurrentPoint');    % button down detected
    finalRect = rbbox;                   % return figure units
    point2 = get(gca,'CurrentPoint');    % button up detected
    point1 = point1(1,1:2);              % extract x and y
    point2 = point2(1,1:2);
    p1 = min(point1,point2);             % calculate locations
    offset = abs(point1-point2);         % and dimensions
%     x = [p1(1) p1(1)+offset(1) p1(1)+offset(1) p1(1) p1(1)];
%     y = [p1(2) p1(2) p1(2)+offset(2) p1(2)+offset(2) p1(2)];
%     hold on
%     axis manual
%     plot(x,y,'r','linewidth',2)          % draw box around selected region
%     hold off;
    
    xmin = floor(p1(1));
    xmax = floor(p1(1) + offset(1));
    ymin = p1(2);
    ymax = p1(2) + offset(2);
%-----------------------------------------------------------------------   
  function[L] = TrimData(Data,L,xmin,xmax,ymin,ymax,flag),
  if xmin < 1 | xmin > length(Data),
      xmin = 1;
  end;
  if xmax > length(Data) | xmax < 1,
      xmax = length(Data);
  end;
  if flag == 1, %selection
      G  =  (Data > ymin & Data < ymax);
      G([1:xmin xmax:length(Data)]) = 0; 
      L = L | G;
  elseif flag == 0, %deselect
      G = ~(Data > ymin & Data < ymax);
      G([1:xmin xmax:length(Data)]) = 1; 
      L = L & G;
  end;
%-----------------------------------------------------------------------